Bayesian Spatio-Temporal modelling of crime against women in Uttar Pradesh (Code)

Abstract

Include your Abstract here. This should be an unstructured abstract of 200 words, summarising the main contents of the scientific report.

Overview

This document provides the code used to carry out the analysis of the project on Bayesian modelling of crime against women in India, for reproducibility purposes. The code is arranged in the following sections, which can be navigated using “Outline” in RMarkdown:

1. Descriptive analysis

2. Modelling: this section includes code used to fit three models: spatial, spatio-temporal, and spatio-temporal with Interaction type I, each fitted separately for rape and dowry deaths, as represented in the subtitles. It also includes the code for visualising the results of the models.

3. Results visualisation: the code used to produce additional figures and compile them into the visualisation file that was used in the report.

# loading all required libraries
library(dplyr)        
library(sf)           
library(spdep)        
library(tidyr)        
library(INLA)         
library(ggplot2)      
library(viridis)      
library(patchwork)    
library(knitr)
library(kableExtra)
library(scales)
library(sp)
library(grid)



# Setting up the workspace
rm(list=ls())
#setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
load("CrimeUttarPradesh.RData")
carto_up <- st_as_sf(carto_up)
ls()

## [1] "carto_india" "carto_up"    "data"

Descriptive analysis

Temporal trends of SMR & visualisation

# Calculate crude rates per 100,000 women - to be used in dataframe below for Supplementary Figure 1
data$rape_rate_per100k <- data$rape / data$pop * 100000
data$dowry_rate_per100k <- data$dowry/ data$pop * 100000

# Create data aggregated across time (i.e. collapsed across space)

year_agg <- data %>% group_by(year) %>% summarise(rape_obs = sum(rape),
                                                  rape_exp = sum(e_rape),
                                                  dowry_obs = sum(dowry),
                                                  dowry_exp = sum(e_dowry),
                                                  total_pop = sum(pop),
                                                  SMR_rape = sum(rape)/sum(e_rape),
                                                  SMR_dowry = sum(dowry)/sum(e_dowry),
                                                  rape_crude_rate = rape_obs / total_pop * 100000,
                                                  dowry_crude_rate = dowry_obs / total_pop * 100000,
                                                  max_rape_region = dist[which.max(rape_rate_per100k)],
                                                  max_rape_rate = max(rape_rate_per100k),
                                                  min_rape_region = dist[which.min(rape_rate_per100k)],
                                                  min_rape_rate = min(rape_rate_per100k),
                                                  max_dowry_region = dist[which.max(dowry_rate_per100k)],
                                                  max_dowry_rate = max(dowry_rate_per100k),
                                                  min_dowry_region = dist[which.min(dowry_rate_per100k)],
                                                  min_dowry_rate = min(dowry_rate_per100k))

# Visualise temporal trends of SMRs

usmoothed_map <- ggplot(year_agg, aes(x = year)) +
  geom_line(aes(y=SMR_rape, color = "Rape")) + 
  geom_line(aes(y=SMR_dowry, color = "Dowry")) + theme_minimal() +
  labs(title = "Temporal evolution of SMRs of Rape and Dowry Crime", 
       x = "Year", y = "Crude Rate", color = "Crime Type")

usmoothed_map

# Visualise summary stats - regions with max & min crude rates each year

year_agg_vis <- year_agg %>% mutate(
  'Year' = year,
  'Rape' = round(rape_crude_rate, 1),
  'Dowry' = round(dowry_crude_rate, 1),
  'Higest Rape' = paste0(max_rape_region, " (", round(max_rape_rate, 1), ")"),
  'Lowest Rape' = paste0(min_rape_region, " (", round(min_rape_rate, 1), ")"),
  'Highest Dowry'= paste0(max_dowry_region, " (", round(max_dowry_rate, 1), ")"),
  'Lowest Dowry' = paste0(min_rape_region, " (", round(min_dowry_rate, 1), ")") ) %>% 
  select( 'Year', 'Rape', 'Dowry', 'Higest Rape', 'Lowest Rape', 'Highest Dowry', 'Lowest Dowry')

reg_table <- knitr::kable(year_agg_vis, caption = "Yearly summary of crude crime rates (per 100,000 women)") %>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "center")
reg_table

Yearly summary of crude crime rates (per 100,000 women)

Year	Rape	Dowry	Higest Rape	Lowest Rape	Highest Dowry	Lowest Dowry
2001	5.4	6.1	Pilibhit (18.4)	Sant Ravidas Nagar Bhadohi (0.3)	Kheri (12.9)	Sant Ravidas Nagar Bhadohi (1.5)
2002	3.8	5.1	Pilibhit (11.3)	Sant Ravidas Nagar Bhadohi (0)	Etah (14.7)	Sant Ravidas Nagar Bhadohi (1.2)
2003	2.4	3.4	Hathras (7.2)	Sant Kabir Nagar (0)	Mainpuri (8.7)	Sant Kabir Nagar (0.8)
2004	3.5	4.3	Sitapur (8.6)	Ghazipur (0.6)	Firozabad (10)	Ghazipur (0.8)
2005	3.0	3.8	Aligarh (8.4)	Ghazipur (0.3)	Firozabad (9.7)	Ghazipur (0.4)
2006	3.1	4.3	Gautam Buddha Nagar (7.7)	Sant Ravidas Nagar Bhadohi (0.6)	Kanpur Dehat (13.4)	Sant Ravidas Nagar Bhadohi (1.6)
2007	3.8	4.8	Chitrakoot (12.1)	Deoria (0.3)	Mainpuri (12.2)	Deoria (1.5)
2008	4.2	5.1	Chitrakoot (11.2)	Sant Ravidas Nagar Bhadohi (0.6)	Aligarh (11.1)	Sant Ravidas Nagar Bhadohi (1.1)
2009	3.9	4.9	Chitrakoot (9.9)	Ghazipur (0.4)	Aligarh (10.2)	Ghazipur (1.9)
2010	3.3	4.8	Rampur (12)	Sant Ravidas Nagar Bhadohi (0.3)	Aligarh (11.3)	Sant Ravidas Nagar Bhadohi (1.1)
2011	4.3	4.9	Aligarh (10.2)	Sant Ravidas Nagar Bhadohi (0.5)	Etawah (11.2)	Sant Ravidas Nagar Bhadohi (1.3)
2012	4.1	4.6	Chitrakoot (14.2)	Mirzapur (0.7)	Aligarh (11)	Mirzapur (1)
2013	6.2	4.7	Hamirpur (19.3)	Mirzapur (0.8)	Firozabad (10.7)	Mirzapur (1.2)
2014	6.8	4.9	Banda (17.3)	Shrawasti (2.1)	Etah (11.8)	Shrawasti (1.3)

Temporal trends visualisation per region

Produced for supplementary figure 1

# Define color map
color_map <- c(
  "Dowry: State Average" = "#8E44AD",
  "Rape: State Average" = "#21918C",
  "Individual Districts" = "lightgray"
)

# Dowry crime 
p_dowry <- ggplot() + 
  geom_line(data=data, aes(x=year, y=dowry_rate_per100k, group=ID_area), 
            color = "lightgray", alpha=0.7, show.legend = FALSE) +
  geom_line(data=year_agg, aes(x=year, y=dowry_crude_rate, color = "Dowry: State Average"), size = 1.5) + 
  scale_color_manual(name = "Legend", values = color_map) +
  theme(legend.position = "right") +
  labs( y = "Crude Rate of Dowry Crime", x = "Year") + 
  theme_minimal()

# Rape 

p_rape <- ggplot() + 
  geom_line(data=data, aes(x=year, y=rape_rate_per100k, group=ID_area), 
            color = "lightgray", alpha=0.7, show.legend = FALSE) +
  geom_line(data=year_agg, aes(x=year, y=rape_crude_rate, color = "Rape: State Average"), size = 1.5) + 
  scale_color_manual(name = "Legend", values = color_map) +
  theme(legend.position = "right") +
  labs(y = "Crude Rate of Rape", x = "Year") + 
  theme_minimal()

# Compule into  single plot 
region_temp_smrs <- ((p_dowry | p_rape) +
  plot_layout(guides = "collect") +
  plot_annotation(title = "Crude Rates of Rape and Dowry Crime in Uttar Pradesh (2001–2014)"))

region_temp_smrs

Spatial visualisation of SMRs

# Aggregate data across space & calculate required summary stats

dist_rates <- data %>%
  group_by(ID_area) %>%
  summarise(
    rape_obs = sum(rape),
    rape_exp = sum(e_rape),
    dowry_obs = sum(dowry),
    dowry_exp = sum(e_dowry),
    SMR_rape = sum(rape)/sum(e_rape),
    SMR_dowry = sum(dowry)/sum(e_dowry))

# Join with SP object 
map_rates <- left_join(carto_up, dist_rates, by = c("ID_area" = "ID_area"))

# Convert to long format for easier plotting
map_rates_long_SMR <- map_rates %>% 
  select(ID_area, geometry, SMR_rape, SMR_dowry) %>%
  pivot_longer(cols = c(SMR_rape, SMR_dowry), names_to = "CrimeType", values_to = "SMR") %>%
  mutate(CrimeType = recode(CrimeType, rape_rate = "Rape",dowry_rate = "Dowry")
  )

# Plot SMRs 
unsmoothed_smr_map <- ggplot(map_rates_long_SMR) +
  geom_sf(aes(fill = SMR), color = NA) +
  scale_fill_viridis_c(
    name = "SMR",
    labels = label_comma(accuracy = 1)) +
  facet_wrap(~ CrimeType) + # set ncol=1 argument if want top and bottom
  labs(title = "SMRs by District (2001–2014)") + theme_bw()

unsmoothed_smr_map

Modelling:

# 1 - Check the coordinates & min distances between polygons:

#st_coordinates(st_centroid(carto_up)) # these appear to be latitudes & longitudes
coord <- st_coordinates(st_centroid(carto_up))
dists <- spDists(coord)
#min(dists[dists > 0]) # min distance is 0.309

# 2- define neighbours list & convert to adjacency matrix

carto_up$reid <- carto_up$ID_area
carto_up_nb <- poly2nb(carto_up, snap=0.001, queen=TRUE) #  snap=0.001 seems the most accurate geographically
summary(carto_up_nb)

## Neighbour list object:
## Number of regions: 70 
## Number of nonzero links: 344 
## Percentage nonzero weights: 7.020408 
## Average number of links: 4.914286 
## Link number distribution:
## 
##  1  2  3  4  5  6  7  8  9 
##  2  1 12 14 16 14  6  3  2 
## 2 least connected regions:
## Uttar Pradesh_Lalitpur Uttar Pradesh_Saharanpur with 1 link
## 2 most connected regions:
## Uttar Pradesh_Budaun Uttar Pradesh_Fatehpur with 9 links

nb2INLA("map_adj",carto_up_nb) # create object with the location of the graph
adj = inla.read.graph(filename="map_adj") # store graph

Spatial model

Dowry

# create dataset collaped across time
up_agg <- data %>% group_by(ID_area) %>% # aggregate over areas 
  summarise(O = sum(dowry), E = sum(e_dowry))  

# fit hierarchical poisson log-linear model
ID <- seq(1,70)

# define formula
f_BYM2 <- O ~ f(ID, model="bym2", graph=adj,
                      hyper=list(
                        prec=list(prior = "pc.prec", param = c(0.5 / 0.31, 0.01)),
                        phi=list(prior="pc", param=c(0.5,2/3))))
# fit model
s_dowry_model = inla(formula=f_BYM2, family="poisson", data=up_agg, E=E, 
               control.predictor = list(compute=TRUE), 
               control.compute = list(waic=TRUE))


# obtain posterior RRs

RR_s <- c()
for(i in 1:70) {RR_s[i] <- inla.emarginal(function(x) exp(x), s_dowry_model$marginals.random$ID[[i]])}

# obtain posterior probabilities

RR_s_marg <- s_dowry_model$marginals.random$ID[1:70]
PP_s <- lapply(RR_s_marg, function(x) {1-inla.pmarginal(0,x)})

# Combine RRs & PPs into a dataframe

s_RR_PP <- data.frame(resRR=RR_s, PP=unlist(PP_s), SP_ID=up_agg[,1])

# Plot posterior RRs & PPs:

# 1 - Categorise variables 
breaks_rr = c(min(s_RR_PP$resRR), 0.6, 0.8, 1, 1.2, 1.4, 1.6, 1.8, 2, max(s_RR_PP$resRR))
s_RR_PP$resRR_cat = cut(s_RR_PP$resRR, breaks=breaks_rr, include.lowest = TRUE)
breaks_pp = c(0, 0.2, 0.8, 1.00)
s_RR_PP$PP_cat = cut(s_RR_PP$PP, breaks_pp, include.lowest = TRUE)

# 2 - join with SP
map_s_RR_PP=left_join(carto_up, s_RR_PP, by=c("ID_area"="ID_area"))

# 3 - plot RRs
dowry_RR_s = ggplot() + geom_sf(data=map_s_RR_PP) + aes(fill=resRR_cat) + theme_bw() +
  scale_fill_brewer(palette = "PuOr") + guides(fill=guide_legend(title="RR")) + ggtitle("RR Spatial model") +
    theme(text = element_text(size=15), 
                  axis.text.x = element_blank(), 
                  axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold"))

# 4 - plot PPs
dowry_PP_s = ggplot() + geom_sf(data=map_s_RR_PP) + aes(fill=PP_cat) + theme_bw() +
    scale_fill_viridis(
    option = "plasma", name="PP",
    discrete = T,
    direction = -1,
    guide = guide_legend(
      title.position = 'top',
      reverse = T
    )) +  ggtitle("PP Spatial model") + theme(text = element_text(size=15), 
                  axis.text.x = element_blank(), 
                  axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold")) 

Rape

# create dataset collaped across time
up_agg_rape <- data %>% group_by(ID_area) %>% # aggregate over areas 
  summarise(O = sum(rape), E = sum(e_rape))

# fit hierarchical poisson log-linear model
ID <- seq(1,70)
f_BYM2 <- O ~ f(ID, model="bym2", graph=adj,
                      hyper=list(
                        prec=list(prior = "pc.prec", param = c(0.5 / 0.31, 0.01)),
                        phi=list(prior="pc", param=c(0.5,2/3))))

s_rape_model = inla(formula=f_BYM2, family="poisson", data=up_agg_rape, E=E, 
               control.predictor = list(compute=TRUE), 
               control.compute = list(waic=TRUE))


# obtain posterior RRs

RR_s <- c()
for(i in 1:70) {RR_s[i] <- inla.emarginal(function(x) exp(x), s_rape_model$marginals.random$ID[[i]])}

# obtain posterior PPs 

RR_s_marg <- s_rape_model$marginals.random$ID[1:70]
PP_s <- lapply(RR_s_marg, function(x) {1-inla.pmarginal(0,x)})

# combine into a dataframe
s_RR_PP <- data.frame(resRR=RR_s, PP=unlist(PP_s), SP_ID=up_agg[,1])

# Plot - as detailed in the chunk above

breaks_rr = c(min(s_RR_PP$resRR), 0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 1.8, 2, max(s_RR_PP$resRR))
s_RR_PP$resRR_cat = cut(s_RR_PP$resRR, breaks=breaks_rr, include.lowest = TRUE)
breaks_pp = c(0, 0.2, 0.8, 1.00)
s_RR_PP$PP_cat = cut(s_RR_PP$PP, breaks_pp, include.lowest = TRUE)
map_s_RR_PP=left_join(carto_up, s_RR_PP, by=c("ID_area"="ID_area"))

rape_RR_s = ggplot() + geom_sf(data=map_s_RR_PP) + aes(fill=resRR_cat) + theme_bw() +
  scale_fill_brewer(palette = "PuOr") + guides(fill=guide_legend(title="RR")) + ggtitle("RR Spatial model") +
    theme(text = element_text(size=15), 
                  axis.text.x = element_blank(), 
                  axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold"))

rape_PP_s = ggplot() + geom_sf(data=map_s_RR_PP) + aes(fill=PP_cat) + theme_bw() +
    scale_fill_viridis(
    option = "plasma", name="PP",
    discrete = T,
    direction = -1,
    guide = guide_legend(
      title.position = 'top',
      reverse = T
    )) +  ggtitle("PP Spatial model") + theme(text = element_text(size=15), 
                  axis.text.x = element_blank(), 
                  axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold")) 

ST model

Dowry

# Prepare data 
data_st_dowry = left_join(carto_up, data, by = "ID_area")
data_st_dowry = data_st_dowry %>% dplyr::rename(O=dowry, E=e_dowry)

# Define formula & fit the model

formula_st_dowry = O ~ f(ID_area, model="bym2", graph=adj, hyper=list(
                        prec=list(prior = "pc.prec", param = c(0.5 / 0.31, 0.01)),
                        phi=list(prior="pc", param=c(0.5,2/3))
                      )) + 
  f(ID_year, model="rw1", hyper=list(prec=list(prior="pc.prec", param=c(0.5/0.31, 0.01))))

st_dowry_model = inla(formula=formula_st_dowry, family="poisson", data=data_st_dowry, E=E, control.compute=list(waic=TRUE), 
                control.predictor = list(compute=TRUE))

# Derive residual posterior means & CIs of random walk component
# Note: they are not plotted yet as they will be first combined with ST+I resRRs - please see "Model Visualisation"

RR_stRW_RR = c()
RR_stRW_l = c()
RR_stRW_h = c()

for(i in 1:14) {
  RR_stRW_RR[i] = inla.emarginal(function(x) exp(x), st_dowry_model$marginals.random$ID_year[[i]])
  RR_stRW_l[i] = inla.qmarginal(0.025,inla.tmarginal(function(x) exp(x), st_dowry_model$marginals.random$ID_year[[i]]))
  RR_stRW_h[i] = inla.qmarginal(0.975,inla.tmarginal(function(x) exp(x), st_dowry_model$marginals.random$ID_year[[i]]))
}

RR_stRW_dowry = data.frame(RR=RR_stRW_RR, low_CI=RR_stRW_l, up_CI = RR_stRW_h)


# Derive & plot RRs & PPs of the BYM2 component

RR_stBYM = c()
for(i in 1:70) {RR_stBYM[i] = inla.emarginal(function(x) exp(x), st_dowry_model$marginals.random$ID_area[[i]])} 

RR_stBYM_marg = st_dowry_model$marginals.random$ID_area[1:70]
PP_stBYM = lapply(RR_stBYM_marg, function(x) {1-inla.pmarginal(0,x)})

resRR_PP_st = data.frame(RR=RR_stBYM, PP=unlist(PP_stBYM), ID_area=carto_up$ID_area)

# plot (using same process as for spatial plots above)

breaks = c(min(resRR_PP_st$RR), 0.6, 0.8, 1, 1.2, 1.4, 1.6, 1.8, 2, max(resRR_PP_st$RR))
resRR_PP_st$RR_cat = cut(resRR_PP_st$RR, breaks=breaks, include.lowest = TRUE)

breaks_pp = c(0, 0.2, 0.8, 1.00)
resRR_PP_st$PP_cat = cut(resRR_PP_st$PP, breaks_pp, include.lowest = TRUE)
map_st_RR_PP=left_join(carto_up, resRR_PP_st, by=c("ID_area"="ID_area"))

# Plot RRs
dowry_RR_st = ggplot() + geom_sf(data=map_st_RR_PP) + aes(fill=RR_cat) + theme_bw() +
  scale_fill_brewer(palette = "PuOr") + guides(fill=guide_legend(title="RR")) + ggtitle("RR Spatio-temporal model") +
    theme(text = element_text(size=15), 
                  axis.text.x = element_blank(), 
                  axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold"))

# Plot PPs
dowry_PP_st = ggplot() + geom_sf(data=map_st_RR_PP) + aes(fill=PP_cat) + theme_bw() +
    scale_fill_viridis(
    option = "plasma", name="PP",
    discrete = T,
    direction = -1,
    guide = guide_legend(
      title.position = 'top',
      reverse = T
    )) +  ggtitle("PP Spatio-temporal model") + theme(text = element_text(size=15), 
                  axis.text.x = element_blank(), 
                  axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold")) 

dowry_RR_st | dowry_PP_st 

ST model - rape

# Prep data 
data_st_rape = left_join(carto_up, data, by = "ID_area")
data_st_rape = data_st_rape %>% dplyr::rename(O=rape, E=e_rape)

# Define & fit the model

formula_st_rape = O ~ f(ID_area, model="bym2", graph=adj, hyper=list(
                        prec=list(prior = "pc.prec", param = c(0.5 / 0.31, 0.01)),
                        phi=list(prior="pc", param=c(0.5,2/3))
                      )) + 
  f(ID_year, model="rw1", hyper=list(prec=list(prior="pc.prec", param=c(0.5/0.31, 0.01))))

st_rape_model = inla(formula=formula_st_rape, family="poisson", data=data_st_rape, E=E, control.compute=list(waic=TRUE), 
                control.predictor = list(compute=TRUE))

# Repeat the process above for plotting of spatial RRs & PPs and (eventually) of temporal RRs

RR_stRW_RR = c()
RR_stRW_l = c()
RR_stRW_h = c()

for(i in 1:14) {
  RR_stRW_RR[i] = inla.emarginal(function(x) exp(x), st_rape_model$marginals.random$ID_year[[i]])
  RR_stRW_l[i] = inla.qmarginal(0.025,inla.tmarginal(function(x) exp(x), st_rape_model$marginals.random$ID_year[[i]]))
  RR_stRW_h[i] = inla.qmarginal(0.975,inla.tmarginal(function(x) exp(x), st_rape_model$marginals.random$ID_year[[i]]))
}

RR_stRW_rape = data.frame(RR=RR_stRW_RR, low_CI=RR_stRW_l, up_CI = RR_stRW_h)


RR_stBYM = c()
for(i in 1:70) {RR_stBYM[i] = inla.emarginal(function(x) exp(x), st_rape_model$marginals.random$ID_area[[i]])} ##
RR_stBYM_marg = st_rape_model$marginals.random$ID_area[1:70]
PP_stBYM = lapply(RR_stBYM_marg, function(x) {1-inla.pmarginal(0,x)})

resRR_PP_st = data.frame(RR=RR_stBYM, PP=unlist(PP_stBYM), ID_area=carto_up$ID_area)

# plot (using same process as described above)
breaks = c(min(resRR_PP_st$RR), 0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 1.8, 2, max(resRR_PP_st$RR))
resRR_PP_st$RR_cat = cut(resRR_PP_st$RR, breaks=breaks, include.lowest = TRUE)
breaks_pp = c(0, 0.2, 0.8, 1.00)
resRR_PP_st$PP_cat = cut(resRR_PP_st$PP, breaks_pp, include.lowest = TRUE)
map_st_RR_PP=left_join(carto_up, resRR_PP_st, by=c("ID_area"="ID_area"))

# RR plot
rape_RR_st = ggplot() + geom_sf(data=map_st_RR_PP) + aes(fill=RR_cat) + theme_bw() +
  scale_fill_brewer(palette = "PuOr") + guides(fill=guide_legend(title="RR")) + ggtitle("RR Spatio-temporal model") +
    theme(text = element_text(size=15), 
                  axis.text.x = element_blank(), 
                  axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold"))

# PP plot

rape_PP_st = ggplot() + geom_sf(data=map_st_RR_PP) + aes(fill=PP_cat) + theme_bw() +
    scale_fill_viridis(
    option = "plasma", name="PP",
    discrete = T,
    direction = -1,
    guide = guide_legend(
      title.position = 'top',
      reverse = T
    )) +  ggtitle("PP Spatio-temporal model") + theme(text = element_text(size=15), 
                  axis.text.x = element_blank(), 
                  axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold")) 
rape_RR_st | rape_PP_st

Spatio-temporal model with Interaction I

Dowry

# Define formula

formula_stInt_dowry = O ~ f(ID_area, model="bym2", graph=adj, hyper=list(
                        prec=list(prior = "pc.prec", param = c(0.5 / 0.31, 0.01)),
                        phi=list(prior="pc", param=c(0.5,2/3))
                      )) + 
  f(ID_year, model="rw1", hyper=list(prec=list(prior="pc.prec", param=c(0.5/0.31, 0.01)))) + 
  f(ID_area_year, model="iid", hyper=list(prec=list(prior="pc.prec", param = c(0.5/0.31, 0.01))))

# Fit model
st_int_model_dowry = inla(formula = formula_stInt_dowry, family="poisson", data=data_st_dowry, E=E, 
                    control.compute=list(dic=TRUE, waic=TRUE, config = TRUE, return.marginals.predictor = TRUE), 
                    control.predictor = list(compute=TRUE))

# extract temporal RRs & CIs to be usied for plotting later

RR_stInt_RW_RR = c()
RR_stInt_RW_l = c()
RR_stInt_RW_h = c()

for(i in 1:14) {
  RR_stInt_RW_RR[i] = inla.emarginal(function(x) exp(x), st_int_model_dowry$marginals.random$ID_year[[i]])
  RR_stInt_RW_l[i] = inla.qmarginal(0.025,inla.tmarginal(function(x) exp(x), st_int_model_dowry$marginals.random$ID_year[[i]]))
  RR_stInt_RW_h[i] = inla.qmarginal(0.975,inla.tmarginal(function(x) exp(x), st_int_model_dowry$marginals.random$ID_year[[i]]))
}

RR_stInt_RW_dowry = data.frame(RR=RR_stInt_RW_RR, low_CI=RR_stInt_RW_l, up_CI = RR_stInt_RW_h)

# extract spatial RRs

RR_stInt_BYM = c()
for(i in 1:70) {RR_stInt_BYM[i] = inla.emarginal(function(x) exp(x), 
                                                  st_int_model_dowry$marginals.random$ID_area[[i]])} 
RR_stInt_BYM_marg = st_int_model_dowry$marginals.random$ID_area[1:70]
PP_stInt_BYM = lapply(RR_stInt_BYM_marg, function(x) {1-inla.pmarginal(0,x)})

# combine 

resRR_PP_stInt = data.frame(RR=RR_stInt_BYM, PP=unlist(PP_stInt_BYM), ID_area=carto_up$ID_area)

# plot (using the process defined above)

breaks = c(min(resRR_PP_stInt$RR), 0.6, 0.8, 1, 1.2, 1.4, 1.6, 1.8, 2, max(resRR_PP_stInt$RR))
resRR_PP_stInt$RR_cat = cut(resRR_PP_stInt$RR, breaks=breaks, include.lowest = TRUE)
breaks_pp = c(0, 0.2, 0.8, 1.00)
resRR_PP_stInt$PP_cat = cut(resRR_PP_stInt$PP, breaks_pp, include.lowest = TRUE)
map_stInt_RR_PP=left_join(carto_up, resRR_PP_stInt, by=c("ID_area"="ID_area"))

#  RR plot
dowry_RR_stInt = ggplot() + geom_sf(data=map_stInt_RR_PP) + aes(fill=RR_cat) + theme_bw() +
  scale_fill_brewer(palette = "PuOr") + guides(fill=guide_legend(title="RR")) + ggtitle("RR Spatio-temporal model with Interaction I") +
    theme(text = element_text(size=15), 
                  axis.text.x = element_blank(), 
                  axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold"))

# PP plot
dowry_PP_stInt = ggplot() + geom_sf(data=map_stInt_RR_PP) + aes(fill=PP_cat) + theme_bw() +
    scale_fill_viridis(
    option = "plasma", name="PP",
    discrete = T,
    direction = -1,
    guide = guide_legend(
      title.position = 'top',
      reverse = T
    )) +  ggtitle("PP Spatio-temporal model with Interaction I") + theme(text = element_text(size=15), 
                  axis.text.x = element_blank(), 
                  axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold")) 

dowry_RR_stInt | dowry_PP_stInt

Spatio-temporal model with Interaction I - rape

formula_stInt_rape = O ~ f(ID_area, model="bym2", graph=adj, hyper=list(
                        prec=list(prior = "pc.prec", param = c(0.5 / 0.31, 0.01)),
                        phi=list(prior="pc", param=c(0.5,2/3))
                      )) + 
  f(ID_year, model="rw1", hyper=list(prec=list(prior="pc.prec", param=c(0.5/0.31, 0.01)))) + 
  f(ID_area_year, model="iid", hyper=list(prec=list(prior="pc.prec", param = c(0.5/0.31, 0.01))))

st_int_model_rape = inla(formula = formula_stInt_rape, family="poisson", data=data_st_rape, E=E, 
                    control.compute=list(dic=TRUE, waic=TRUE))

# temporal RRs & CIs

RR_stInt_RW_RR = c()
RR_stInt_RW_l = c()
RR_stInt_RW_h = c()

for(i in 1:14) {
  RR_stInt_RW_RR[i] = inla.emarginal(function(x) exp(x), st_int_model_rape$marginals.random$ID_year[[i]])
  RR_stInt_RW_l[i] = inla.qmarginal(0.025,inla.tmarginal(function(x) exp(x), st_int_model_rape$marginals.random$ID_year[[i]]))
  RR_stInt_RW_h[i] = inla.qmarginal(0.975,inla.tmarginal(function(x) exp(x), st_int_model_rape$marginals.random$ID_year[[i]]))
}

RR_stInt_RW_rape = data.frame(RR=RR_stInt_RW_RR, low_CI=RR_stInt_RW_l, up_CI = RR_stInt_RW_h)

# spatial RRs

RR_stInt_BYM = c()
for(i in 1:70) {RR_stInt_BYM[i] = inla.emarginal(function(x) exp(x), 
                                                  st_int_model_rape$marginals.random$ID_area[[i]])} 
RR_stInt_BYM_marg = st_int_model_rape$marginals.random$ID_area[1:70]
PP_stInt_BYM = lapply(RR_stInt_BYM_marg, function(x) {1-inla.pmarginal(0,x)})

# paste

resRR_PP_stInt = data.frame(RR=RR_stInt_BYM, PP=unlist(PP_stInt_BYM), ID_area=carto_up$ID_area)


breaks = c(min(resRR_PP_stInt$RR), 0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 1.8, 2, max(resRR_PP_stInt$RR))
resRR_PP_stInt$RR_cat = cut(resRR_PP_stInt$RR, breaks=breaks, include.lowest = TRUE)
breaks_pp = c(0, 0.2, 0.8, 1.00)
resRR_PP_stInt$PP_cat = cut(resRR_PP_stInt$PP, breaks_pp, include.lowest = TRUE)
map_stInt_RR_PP=left_join(carto_up, resRR_PP_stInt, by=c("ID_area"="ID_area"))

c = ggplot() + geom_sf(data=map_stInt_RR_PP) + aes(fill=RR_cat) + theme_bw() +
  scale_fill_brewer(palette = "PuOr") + guides(fill=guide_legend(title="RR")) + ggtitle("RR Spatio-temporal model with Interaction I") +
    theme(text = element_text(size=15), 
                  axis.text.x = element_blank(), 
                  axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold"))

rape_RR_stInt = ggplot() + geom_sf(data=map_stInt_RR_PP) + aes(fill=RR_cat) + theme_bw() +
    scale_fill_viridis(
    option = "plasma", name="PP",
    discrete = T,
    direction = -1,
    guide = guide_legend(
      title.position = 'top',
      reverse = T
    )) +  ggtitle("PP Spatio-temporal model with Interaction I") + theme(text = element_text(size=15), 
                  axis.text.x = element_blank(), 
                  axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold")) 

rape_PP_stInt = ggplot() + geom_sf(data=map_stInt_RR_PP) + aes(fill=PP_cat) + theme_bw() +
    scale_fill_viridis(
    option = "plasma", name="PP",
    discrete = T,
    direction = -1,
    guide = guide_legend(
      title.position = 'top',
      reverse = T
    )) +  ggtitle("PP Spatio-temporal model with Interaction I") + theme(text = element_text(size=15), 
                  axis.text.x = element_blank(), 
                  axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold")) 

rape_RR_stInt | rape_PP_stInt

Results visualisation

Additional figures

Temporal resRRs & CIs (rape & dowry)

# Spatio-temporal model - prepare temporal residual RRs 
RR_stRW_dowry$year <- seq(2001, 2014)
RR_stRW_rape$year   <- seq(2001, 2014)
RR_stRW_dowry$Crime <- "Dowry Deaths"
RR_stRW_rape$Crime   <- "Rape"

RR_stRW_both <- rbind(RR_stRW_dowry, RR_stRW_rape) # Combine 

# Spatio-temporal + Int model - prepare temporal residual RRs 

RR_stInt_RW_dowry$year <- seq(2001, 2014)
RR_stInt_RW_rape$year   <- seq(2001, 2014)
RR_stInt_RW_dowry$Crime <- "Dowry Deaths"
RR_stInt_RW_rape$Crime   <- "Rape"

RR_stInt_RW_both <- rbind(RR_stInt_RW_dowry, RR_stInt_RW_rape) # Combine

# Plot - ST model, both crimes

temp_st_both <- ggplot(RR_stRW_both, aes(x = year, y = RR, colour = Crime, fill = Crime)) +
  geom_line() +
  geom_ribbon(aes(ymin = low_CI, ymax = up_CI), alpha = 0.2, colour = NA) +
  ggtitle("Spatio-temporal model") +
  theme_bw() +
  labs(x = "Year")

# Plot - ST+I model, both crimes

temp_stInt_both = ggplot(RR_stInt_RW_both, aes(x = year, y = RR, colour = Crime, fill = Crime)) +
  geom_line() +
  geom_ribbon(aes(ymin = low_CI, ymax = up_CI), alpha = 0.2, colour = NA) +
  ggtitle("Spatio-temporal model with Interaction I") +
  theme_bw() +
  labs(x = "Year")

# Combine into one figure
temp_st_stInt <- temp_st_both | temp_stInt_both 
temp_st_stInt

Interactions - spatial plot

For both models

# Prepare data 

data_st = left_join(carto_up, data, by = "ID_area")

# Extract interactions 

data_st$int_dowry = st_int_model_dowry$summary.random$ID_area_year$mean
data_st$int_dowry_cat = cut(data_st$int_dowry,  breaks=c(-1,-0.05, 
                  -0.01, 0.01, 0.05, 1),include.lowest = T)

data_st$int_rape = st_int_model_rape$summary.random$ID_area_year$mean
data_st$int_rape_cat = cut(data_st$int_rape,  breaks=c(-1,-0.05, 
                  -0.01, 0.01, 0.05, 1),include.lowest = T)

# Plot interactions for dowry model
dowry_int_map =  ggplot() + geom_sf(data = data_st, aes(fill=int_dowry_cat)) + theme_bw() +
  scale_fill_brewer(palette = "PuOr") +
  guides(fill=guide_legend(title=NULL)) + 
  theme(text = element_text(size=20),
        axis.text.x = element_blank(), axis.text.y = element_blank()) +
facet_wrap(~ year, ncol = 5, labeller=labeller(ID_year=c("1"="2001","2"="2002","3"="2003","4"="2004","5"="2005",
                            "6"="2006","7"="2007","8"="2008","9"="2009","10"="2010",
                            "11"="2011","12"="2012","13"="2013","14"="2014"))) +
labs(title="Dowry")

# Plot interactions for rape model
rape_int_map =  ggplot() + geom_sf(data = data_st, aes(fill=int_rape_cat)) + theme_bw() +
  scale_fill_brewer(palette = "PuOr") +
  guides(fill=guide_legend(title=NULL)) + 
  theme(text = element_text(size=20),
        axis.text.x = element_blank(), axis.text.y = element_blank()) +
facet_wrap(~ year, ncol = 5, labeller=labeller(ID_year=c("1"="2001","2"="2002","3"="2003","4"="2004","5"="2005",
                            "6"="2006","7"="2007","8"="2008","9"="2009","10"="2010",
                            "11"="2011","12"="2012","13"="2013","14"="2014"))) +
labs(title="Rape")

# Combine
int_map = rape_int_map / dowry_int_map 
int_map

Total posterior RRs from ST+I model

Dowry

# Exract posterior RRs

data_st$RR_total_emarg_dowry <- mapply(function(area_idx, year_idx, area_year_idx) {
  rr_area <- inla.emarginal(function(x) exp(x), st_int_model_dowry$marginals.random$ID_area[[area_idx]])
  rr_year <- inla.emarginal(function(x) exp(x), st_int_model_dowry$marginals.random$ID_year[[year_idx]])
  rr_interaction <- inla.emarginal(function(x) exp(x), st_int_model_dowry$marginals.random$ID_area_year[[area_year_idx]])
  
  rr_area * rr_year * rr_interaction
},
area_idx = data_st$ID_area,
year_idx = data_st$ID_year,
area_year_idx = data_st$ID_area_year)

# categorise 
breaks_total = c(min(data_st$RR_total_emarg_dowry), 0.5, 1, 1.5, 2, 2.5, max(data_st$RR_total_emarg_dowry))
data_st$RR_total_cat_dowry = cut(data_st$RR_total_emarg_dowry,  breaks=breaks_total,include.lowest = T)

# Plot for dowry
dowry_stInt_postRR = ggplot() + geom_sf(data = data_st, aes(fill = RR_total_cat_dowry)) + theme_bw() + scale_fill_brewer(palette = "PuOr") +
  guides(fill=guide_legend(title=NULL)) +
  theme(text = element_text(size=20),
        axis.text.x = element_blank(), axis.text.y = element_blank()) +
facet_wrap(~ year, ncol = 5, labeller=labeller(ID_year=c("1"="2001","2"="2002","3"="2003","4"="2004","5"="2005",
                            "6"="2006","7"="2007","8"="2008","9"="2009","10"="2010",
                            "11"="2011","12"="2012","13"="2013","14"="2014"))) +
labs(title = "Dowry Crime")

Rape

data_st$RR_total_emarg_rape <- mapply(function(area_idx, year_idx, area_year_idx) {
  rr_area <- inla.emarginal(function(x) exp(x), st_int_model_rape$marginals.random$ID_area[[area_idx]])
  rr_year <- inla.emarginal(function(x) exp(x), st_int_model_rape$marginals.random$ID_year[[year_idx]])
  rr_interaction <- inla.emarginal(function(x) exp(x), st_int_model_rape$marginals.random$ID_area_year[[area_year_idx]])
  
  rr_area * rr_year * rr_interaction
},
area_idx = data_st$ID_area,
year_idx = data_st$ID_year,
area_year_idx = data_st$ID_area_year)


breaks_total = c(min(data_st$RR_total_emarg_rape), 0.5, 1, 1.5, 2, 2.5, max(data_st$RR_total_emarg_rape))
data_st$RR_total_cat_rape = cut(data_st$RR_total_emarg_rape,  breaks=breaks_total,include.lowest = T)

rape_stInt_postRR = ggplot() + geom_sf(data = data_st, aes(fill = RR_total_cat_rape)) + theme_bw() + scale_fill_brewer(palette = "PuOr") +
  guides(fill=guide_legend(title=NULL)) +
  theme(text = element_text(size=20),
        axis.text.x = element_blank(), axis.text.y = element_blank()) +
facet_wrap(~ year, ncol = 5, labeller=labeller(ID_year=c("1"="2001","2"="2002","3"="2003","4"="2004","5"="2005",
                            "6"="2006","7"="2007","8"="2008","9"="2009","10"="2010",
                            "11"="2011","12"="2012","13"="2013","14"="2014"))) +
labs(title = "Rape")

# Combine posterior RRs & plot 
meanRR_post  = rape_stInt_postRR / dowry_stInt_postRR + plot_annotation(title="Posterior Mean Relative Risks", theme=theme(plot.title=element_text(size=24)))
meanRR_post

Combine remaining figures for visualisation (figures not displayed for considerations of file length)

# descriptive analysis 
unsmoothed <- usmoothed_map | unsmoothed_smr_map
region_temp_smrs <- (p_dowry | p_rape) +
  plot_layout(guides = "collect") +
  plot_annotation(title = "Crude Rates of Rape and Dowry Crime in Uttar Pradesh (2001–2014)")

# supplementary figures of all models' smoothed spatial RR maps 
rape_all_m <- (rape_RR_s | rape_PP_s) / (rape_RR_st | rape_PP_st) / (rape_RR_stInt | rape_PP_stInt)
dowry_all_m <- (dowry_RR_s | dowry_PP_s) / (dowry_RR_st | dowry_PP_st) / (dowry_RR_stInt | dowry_PP_stInt)

# ST+I RR & PP plots for rape and dowry crime (model of interest discussed in report)

rape_title <- wrap_elements(full = textGrob("Rape", rot = 90, gp = gpar(fontsize = 15))) # Make titles
dowry_title <- wrap_elements(full = textGrob("Dowry Crime", rot = 90, gp = gpar(fontsize = 15)))

rape_row <- rape_title + (rape_RR_stInt + rape_PP_stInt + plot_layout(ncol = 2)) +plot_layout(widths = c(0.1, 1)) # Compile rows
dowry_row <- dowry_title + (dowry_RR_stInt + dowry_PP_stInt + plot_layout(ncol = 2)) +plot_layout(widths = c(0.1, 1))

stInt_m <- rape_row / dowry_row # compile figures

Tables

WAICs

WAIC_df = data.frame(model = c("Spatio-temporal", "Spatio-Temporal + Int"), 
                       WAIC_dowry = round(c(st_dowry_model$waic$waic, st_int_model_dowry$waic$waic)),
                        WAIC_rape = round(c(st_rape_model$waic$waic, st_int_model_rape$waic$waic)))
colnames(WAIC_df) = c("Model", "WAIC (Dowry)", "WAIC (Rape)")                    
waic_table = knitr::kable(WAIC_df, caption = "Table 1: WAIC for Different Models and Crimes") %>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "center")
waic_table

Table 1: WAIC for Different Models and Crimes

Model	WAIC (Dowry)	WAIC (Rape)
Spatio-temporal	6654	7781
Spatio-Temporal + Int	6318	6240

extract all hyperparameters

# Combine hyperparameters into one table

rape_s_hyper = s_rape_model$summary.hyperpar
dowry_s_hyper = s_dowry_model$summary.hyperpar
rape_st_hyper = st_rape_model$summary.hyperpar
dowry_st_hyper = st_dowry_model$summary.hyperpar
rape_stInt_hyper = st_int_model_rape$summary.hyperpar
dowry_stInt_hyper = st_int_model_dowry$summary.hyperpar

# tidy up into dataframes: spatial
rape_s_df <- as.data.frame(rape_s_hyper) %>%
  mutate(Hyperparameter = rownames(rape_s_hyper),Crime = "Rape", Model = "Spatial")
dowry_s_df <- as.data.frame(dowry_s_hyper) %>%
  mutate(Hyperparameter = rownames(dowry_s_hyper),Crime = "Dowry",Model = "Spatial")

# spatio-temporal
rape_st_df <- as.data.frame(rape_st_hyper) %>%
  mutate(Hyperparameter = rownames(rape_st_hyper),Crime = "Rape",Model = "Spatio-temporal")
dowry_st_df <- as.data.frame(dowry_st_hyper) %>%
  mutate(Hyperparameter = rownames(dowry_st_hyper),Crime = "Dowry",Model = "Spatio-temporal")

# ST + interaction
rape_stInt_df <- as.data.frame(rape_stInt_hyper) %>%
  mutate(Hyperparameter = rownames(rape_stInt_hyper),Crime = "Rape",Model = "Spatio-temporal Interaction")
dowry_stInt_df <- as.data.frame(dowry_stInt_hyper) %>%
  mutate(Hyperparameter = rownames(dowry_stInt_hyper),Crime = "Dowry",Model = "Spatio-temporal Interaction")

# Combine dataframes
hyper_df <- bind_rows(rape_s_df, dowry_s_df,rape_st_df, dowry_st_df,rape_stInt_df, dowry_stInt_df)

# Plot all hyperparameters from all models
hyper_all = ggplot(hyper_df, aes(x = Model, y = mean, color = Crime, group = Crime)) +
  geom_point(position = position_dodge(width = 0.5), size = 3) +
  geom_errorbar(aes(ymin = `0.025quant`, ymax = `0.975quant`),
                width = 0.2, position = position_dodge(width = 0.5)) +
  facet_wrap(~ Hyperparameter, scales = "free_y") +
  labs(title = "Hyperparameter Estimates by Model and Crime Type",
       y = "Estimate", x = "Model") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

hyper_all

ST+I model hyperparameter table

rape_hyper <- as.data.frame(st_int_model_rape$summary.hyperpar)
dowry_hyper <- as.data.frame(st_int_model_dowry$summary.hyperpar)

# 2. Clean / Rename Columns for Clarity
# Typically, you have columns like "mean", "sd", "0.025quant", "0.975quant", etc.
rape_hyper <- rape_hyper %>%
  mutate(Parameter = rownames(rape_hyper),
         CrimeType = "Rape") %>%
  select(Parameter, mean, `0.025quant`, `0.975quant`, CrimeType)

dowry_hyper <- dowry_hyper %>%
  mutate(Parameter = rownames(dowry_hyper),
         CrimeType = "Dowry Deaths") %>%
  select(Parameter, mean, `0.025quant`, `0.975quant`, CrimeType)

# Combine into a single data frame
hyper_combined <- bind_rows(rape_hyper, dowry_hyper)

# 3. Create a Nice Table
# This example shows a side-by-side listing of parameter, mean, and credible intervals
# for each crime type.

hyper_table <- hyper_combined %>%
  mutate(
    # Format the credible interval as a single string
    CrI = paste0("(", round(`0.025quant`, 3), ", ", round(`0.975quant`, 3), ")"),
    Mean = round(mean, 3)
  ) %>%
  select(CrimeType, Parameter, Mean, CrI)

# 1. Pivot to wide format
hyper_wide <- hyper_combined %>%
  pivot_wider(
    id_cols = Parameter,
    names_from = CrimeType,              # Rape vs. Dowry Deaths
    values_from = c(mean, `0.025quant`, `0.975quant`),
    names_glue = "{CrimeType}_{.value}"  # to create columns like Rape_mean, Dowry Deaths_mean, etc.
  )

# 2. Create columns that combine the mean and 95% CrI into one string
hyper_wide <- hyper_wide %>%
  mutate(
    Rape_Est = paste0(
      round(Rape_mean, 3), 
      " (", round(Rape_0.025quant, 3), ", ", round(Rape_0.975quant, 3), ")"
    ),
    Dowry_Est = paste0(
      round(`Dowry Deaths_mean`, 3), 
      " (", round(`Dowry Deaths_0.025quant`, 3), ", ", round(`Dowry Deaths_0.975quant`, 3), ")"
    )
  ) %>%
  # 3. Select and rename columns for clarity
  select(
    Parameter,
    `Rape: Mean (95% CrI)` = Rape_Est,
    `Dowry Deaths: Mean (95% CrI)` = Dowry_Est
  )

# 4. Create a publication-ready table
stInt_hyper = kable(hyper_wide, "html", caption = "Table 2:Posterior Summaries of Hyperparameters") %>%
  kable_styling(full_width = FALSE)

stInt_hyper

Table 2:Posterior Summaries of Hyperparameters

Parameter	Rape: Mean (95% CrI)	Dowry Deaths: Mean (95% CrI)
Precision for ID_area	5.91 (3.877, 8.516)	11.342 (7.555, 16.159)
Phi for ID_area	0.849 (0.571, 0.982)	0.894 (0.675, 0.99)
Precision for ID_year	15.09 (5.964, 30.721)	52.5 (19.614, 110.791)
Precision for ID_area_year	12.55 (10.766, 14.554)	43.712 (34.533, 54.927)

Compile all figures in .RData file

This is done so that figures can be subsequently loaded into the report file and visualised

save(unsmoothed, region_temp_smrs, rape_all_m, dowry_all_m, rape_title, dowry_title, stInt_m, temp_st_stInt, int_map, meanRR_post, waic_table, stInt_hyper, hyper_all, reg_table, file="Figures.RData")